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The critical behavior of long straight rigid rods of length k (fc-mers) on square and triangular- 
lattices at intermediate density has been studied. A nematic phase, characterized by a big domain of 
parallel fc-mers, was found. This ordered phase is separated from the isotropic state by a continuous 
transition occurring at a intermediate density 9 C . Two analytical techniques were combined with 
Monte Carlo simulations to predict the dependence of 9 C on k, being 9 c (k) oc fc _1 . The first involves 
simple geometrical arguments, while the second is based on entropy considerations. Our analysis 
allowed us also to determine the minimum value of k (k m in — 7), which allows the formation of a 
nematic phase on a triangular lattice. 
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I. INTRODUCTION 

The isotropic-nematic (I-N) phase transition has been a topic of active theoretical and 
experimental studies over the past few decades. An early seminal contribution to this 
subject was made by Onsager 1 with his work on the I-N phase transition of infinitely 
thin rods. This theory predicted that excluded volume interactions alone can lead to 
long-range orientational (nematic) order. Later, several papers contributed greatly to the 
understanding of the statistics of rigid rods^^^ 1 ^. Successive works have established 
detailed phase diagrams for several hard-body shapes^&i&^k For the continuum problem, 
there is general agreement that in three dimensions, infinitely thin rods undergo a first- 
order I-N transition, as was pointed out by Onsager-. In two dimensions, the nature of 
the isotropic-nematic (I-N) phase transition depends crucially on the particle interactions 
and a rich variety of behaviors is observecr^ 1 ^. 

A notable feature is that nematic order is only stable for sufficiently large aspect ratios 
while isotropic systems of short rods not show nematic order at all. The long-range ori- 
entational order also disappears in the case of irreversible adsorption (no desorption)^^, 
where the distribution of adsorbed objects is different from that obtained at equilib- 
rium^^. Thus, at high coverage, the equilibrium state corresponds to a nematic phase 
with long-range correlations, whereas the final state generated by irreversible adsorption 
has infinite memory of the process and orientational order is purely local. 

In the case of lattice models, two previous articles^ 1 ^, referred to as papers I and II, 
respectively, were devoted to the study of the I-N phase transition in a system of long 
straight rigid rods of length k (fc-mers) on two-dimensional lattices with discrete allowed 
orientations. The model of a two-dimensional gas of rigid fc-mers is the simplest repre- 
sentation of a strongly adsorbed film of linear molecules in submonolayer or monolayer 



regime. Examples of this kind of systems are monolayer films of n-alkanes adsorbed on 
mono crystalline surfaces of metals, such as Pt(lll)— and Au(lll}^ i 2^. 



Papers I and II were inspired in the excellent paper by Ghosh and Dhar— . In Ref. 23], 
the authors presented strong numerical evidence that a system of square geometry, with 
two allowed orientations, shows nematic order at intermediate densities for k > 7 and 
provided a qualitative description of a second phase transition (from a nematic order to a 
non-nematic state) occurring at a density close to 1. However, the authors were not able 
to determine the critical quantities (critical point and critical exponents) characterizing 
the I-N phase transition occurring in the system. 

In this context, extensive Monte Carlo (MC) simulations were used in paper I to resolve 



the universality class of the first phase transition occurring in the system of Ref. [23]. 
Lattices of various sizes were considered and finite-size scaling theory was utilized. As 
it was evident from the calculation of the critical exponents and the behavior of Binder 
cumulants, the universality class was shown to be that of the 2D Ising model for square 
lattices with two allowed orientations, and the three-state Potts model for triangular 
lattices with three allowed orientations. 

Paper II was a step further, analyzing the I-N phase transition in terms of entropy. For 
this purpose, the configurational entropy of a system of rigid rods deposited on a square 
lattice was calculated by Monte Carlo (MC) simulations and thermodynamic integration 
method^. The numerical data were compared with the corresponding ones obtained from 
a fully aligned system (nematic phase), whose calculation reduces to the one-dimensional 
case2£. The study allowed us 1) to confirm previous results in the literature^, namely, 
the existence of i) a I-N phase transition at intermediate densities for k > 7 and ii) a 
second phase transition from a nematic order to a non-nematic state at high density. In 



the second case, Ref. 19] represents the first numerical evidence existing in the literature 
about this important point; 2) to provide an interpretation on the underlying physics of 
the observed k dependence of the isotropic-nematic phase transition; and 3) to test the 
predictions of the main theoretical models developed to study adsorption with multisite- 
occupancy. 



23|, 



Even though many aspects of the problem have been treated in papers I, II and Ref. 
other points remain still open. Among them, there exists no studies on the dependence 
of the critical density characterizing the I-N phase transition, 9 C , on the size k of the 
rod. In this work we attempt to remedy this situation. For this purpose, extensive MC 
simulations complemented by finite-size scaling techniques and theoretical modeling have 
been applied. We restrict the study to the first phase transition occurring in the system 
(or I-N phase transition at intermediate density). Our analysis allowed us (1) to obtain 
6 C as a function of k for square and triangular lattices, being 9 c (k) oc k~ l ] and (2) to 
determine the minimum value of k (k m i n = 7), which allows the formation of a nematic 
phase on triangular lattices. 

The outline of the paper is as follows. In Sec. II we describe the lattice-gas model, 
the simulation scheme, and we present the behavior of 6 c (k), obtained by using the MC 
method. In Sec. Ill we present the analytical approximations and compare the MC results 
with the theoretical calculations. Finally, the general conclusions are given in Sec. IV. 

II. LATTICE-GAS MODEL AND MONTE CARLO SIMULATION SCHEME 
A. Model and Monte Carlo method 



We address the general case of adsorbates assumed to be linear rigid particles containing 
k identical units (fc-mers), with each one occupying a lattice site. Small adsorbates would 



correspond to the monomer limit (k = 1). The distance between fc-mer units is assumed 
to be equal to the lattice constant; hence exactly k sites are occupied by a k-mer when 
adsorbed. The only interaction between different rods is hard-core exclusion: no site 
can be occupied by more than one k-mer unit. The surface is represented as an array 
of M = L x L adsorptive sites in a square or triangular lattice arrangement, where L 
denotes the linear size of the array. 

The degree of order in the adsorbed phase is calculated for each configuration according 
to the standard method used for the Potts modePi. To this end, we first build a set of 
vectors {f? 1; n 2 , • ■ ■ , n m } with the following properties: (i) each vector is associated to one 
of the m possible orientations (or directions) for a fc-mer on the lattice; (ii) the n^s lie in 
the same plane (or are co-planar) and point radially outward from a given point P which 
is defined as coordinate origin; (in) the angle between two consecutive vectors, n« and 
n i+ i, is equal to 2n/m; and {iv) the magnitude of is equal to the number of fc-mers 
aligned along the i-direction. Note that the n^s have the same directions as the q vectors 
in Ref These directions are not coincident with the allowed directions for the fc-mers 
on the real lattice. Then the order parameter 5 of the system is given by 

S = ^Fp- { (l) 

5 represents a general order parameter measuring the orientation of the fc-mers on a lattice 
with m directions. In the case of square lattices, m = 2 and the angle between n\ and n,2 
is 7r. Accordingly, the order parameter reduces to 5 = \n\ — n 2 \ / {ni + 112), n\ (n 2 ) being 
the number of fc-mers aligned along the horizontal (vertical) direction. This expression 
coincides with the order parameter Q defined in Ref. 2j|. On the other hand, m = 3 for 
triangular lattices and n\, n 2 and 713 form angles of 27r/3 between them. 
When the system is disordered (8 < 6 C ), all orientations are equivalents and 5 is zero. 



As the density is increased above 6 C , the fc-mers align along one direction and 5 is different 
from zero. Thus, 5 appears as a proper order parameter to elucidate the phase transition. 

The problem has been studied by grand canonical MC simulations using a typical 
adsorption-desorption algorithm2L22a22. The procedure is as follows. Once the value of the 
chemical potential /i is set, a linear k-uple of nearest-neighbor sites is chosen at random. 
Then, if the k sites are empty, an attempt is made to deposit a rod with probability 
W = min {1, exp (/x//cbT)}, where ks is the Boltzmann constant and T is the tempera- 
ture; if the k sites are occupied by units belonging to the same fc-mer, an attempt is made 
to desorb this k-mer with probability W = min {1, exp (— /i/^T)}; and otherwise, the 
attempt is rejected. In addition, displacement (diffusional relaxation) of adparticles to 
nearest-neighbor positions, by either jumps along the k-mer axis or reptation by rotation 
around the k-mer end, must be allowed in order to reach equilibrium in a reasonable time. 
A MC step (MCs) is achieved when M fc-uples of sites have been tested to change its oc- 
cupancy state. Typically, the equilibrium state can be well reproduced after discarding 
the first r' = 10 6 MCs. Then, the next r = 2 x 10 6 MCs are used to compute averages. 

In our Monte Carlo simulations, we varied the chemical potential \i and monitored the 
density 9 and the order parameter 5, which can be calculated as simple averages. The 
reduced fourth-order cumulant Ul introduced by Binder— was calculated as: 

U - 1 (0 (2) 

where (• • • ) means the average over the r MC simulation runs. 

Finally, the configurational entropy of the system S, was calculated by using thermody- 
namic integration method JkhMiMi2& m The method in the grand canonical ensemble relies 
upon integration of the chemical potential fi on coverage along a reversible path between 
an arbitrary reference state and the desired state of the system. This calculation also 



requires the knowledge of the total energy U for each obtained coverage. Thus, for a 
system made of N particles on M lattice sites, we have: 

S(N, M, T) — S(N 0; M, T) + M^IhJWM 

i r N 

-- / ^dN' (3) 

In our case U(N, M, T) = and the determination of the entropy in the reference state, 
S(N , M, T), is trivial [S(N ,M,T) = for iV = 0]. Note that the reference state, 
N — > 0, is obtained for \ijksT — > — oo. Then, 

k B k B T J k 

where s(= S/M) is the configurational entropy per site and 6(= k N/M) is the surface 
coverage (or density). 

B. Computational results 

Computational simulations have been developed for a system of straight rigid rods of 
length k (k = 2 — 14) on a lattice. The surface was represented as an array of ad- 
sorptive sites in a square or triangular lattice arrangement. In addition, conventional 
periodic boundary conditions were considered. The effect of finite size was investi- 
gated by examining square lattices with L/k = 5, 10, 15,20 and triangular lattices with 
L/k = 10,15,20,25. 



As it was established in Ref. [23j, the minimum value of k, which allows the formation 
of a nematic phase on a square lattice at intermediate densities, is k m i n = 7. This critical 
quantity has not been calculated yet for other geometries. Then, our first objective is 
to obtain k m i n for triangular lattices. For this purpose, two criteria have been applied: 
1) the comparison between the configurational entropy of the system, obtained by MC 



simulations, and the corresponding to a fully aligned system (nematic phase), whose 
calculation reduces to the ID case^; and 2) the behavior of the nematic order parameter 
function of coverage. 

The results obtained in the first case are shown in Fig. 1. Dotted line and symbols 
represent MC data for triangular lattices, k = 6, 7 and 8 and L/k = 20. The calculation of 
s(6)/kB through eq. (jl]) is straightforward and computationally simple, since the coverage 
dependence of /i/fceT is evaluated following the standard procedure of MC simulation 
described in previous section. Then, fi(9)/kBT is spline-fitted and numerically integrated. 

On the other hand, when the nematic phase is formed, the system is characterized by 
a big domain of parallel fc-mers. The calculation of the entropy of this fully aligned state 
having density 9 reduces to the calculation of a one- dimensional problem 25 
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Results from eq. ([5]) for k = 6, 7 and 8 are shown in Fig. 1 (solid lines). As it can be 
observed, for k < 6, the ID results present a smaller s//cg than the 2D simulation data 
over all the range of 9. For k > 7, there exists a range of coverage for which the difference 
between the ID value and the true 2D value is very small. In other words, for k > 7 and 
intermediate densities, it is more favorable for the rods to align spontaneously, since the 
resulting loss of orientational entropy is by far compensated by the gain of translational 
entropy. This finding is a clear indication that k m i n = 7 for triangular lattices. In addi- 
tion, 9 C can be calculated from the minimum value of 9 for which the near superposition 
of the ID and 2D results occurs. Thus, the technique provides an alternative method 
of determining the critical coverage characterizing the I-N phase transition without any 
special requirement and time consuming computation. However, it is important to em- 



phasize that the calculation of the entropy of the nematic phase from the ID model is an 
approximation (especially at the moderate densities, where the phase is not completely 
aligned). Consequently, a precise determination of 9 C should require an extensive work of 
MC simulation and finite-size scaling techniques. 

As an independent corroboration of the results previously obtained, the inset in Fig. 1 
presents the nematic order parameter 5 as a function of coverage for k = 6, 7 and 8 and 
L/k = 20. The behavior of 5 confirms the existence of nematic order for k > 7. This 
value of k m in coincides with that obtained for square lattices. 

Once k min has been established, we now discuss the behavior of the critical density 
as a function of the size k. In the case of the standard theory of finite-size scaling^^, 
when the phase transition is temperature driven, the technique allows for various efficient 
routes to estimate T c from MC data. One of these methods, which will be used in this 
case, is from the temperature dependence of Ul{T), which is independent of the system 
size for T = T c . In other words, T c is found from the intersection of the curve Ul(T) for 
different values of L, since Ul(T c ) = const. In our study, we modified the conventional 
finite-size scaling analysis by replacing temperature by density^. Under this condition, 
the critical density has been estimated from the plots of the reduced four-order cumulants 
Ul{6) plotted versus 9 for several lattice sizes. As an example, Fig. 2 shows the results 
for rods of size k = 10 on square [Fig. 2(a)] and triangular [Fig. 2(b)] lattices. In the 
cases of the figure, the values obtained were 9 C = 0.502(1) and 9 C = 0.530(1), for square 
and triangular lattices, respectively. The curves of the order parameter, which were used 
to obtain Ul(0), are shown in the insets of the figure. 

The procedure of Fig. 2 was repeated for 7 < k < 14 and the results are collected in 
Fig. 3. The log-log plots show that the critical density follows a power law as 9 c (k) oc k~ l . 



The understanding of this dependence of 9 c (k) on k can be developed by using simple 
geometric arguments. An example follows in order to make this point clear. Fig. 4(a) 
shows a snapshot corresponding to a "small window" of side I, which is embedded in an 
infinite square lattice. Open squares correspond to empty sites and black bars represent 
adsorbed fc-mers (in the case of the figure, k = 4 and 1 = 8). Then, the local coverage or 
fraction of occupied sites belonging to the window is 9 = kN/ (I x I), being N the number 
of fc-mers into the window^. We can now think of a transformation L — > L' from the 
original lattice L to a new lattice L', where each characteristic length of L' is s times the 
corresponding one of L (s = 1,2,3,4, • ■ •). Then, I' = si, k! = sk and 9' = k'N/(l' x I'). 
By following this procedure, Fig. 4(b) was obtained from Fig. 4(a). In this case, s = 3, 
and consequently, I' = 31 = 24 and k' = 3k = 12. As it can be observed, L and V 
represent the "same situation" for two different values of k. 
In general, 

&_ = k'N(l x I) 

9 kN{l' xl'Y U 
and l/k = l'/k'. Then, the relationship between a coverage in L and the corresponding 
one in V results 

- = (7) 
9' k K ) 

By taking the system L' as reference, it is possible to write, 

9 = (k'^k- 1 

oc k~\ (8) 

as it was found from the MC data. A more rigorous justification of the observed k 
dependence of the I-N phase transition will be presented in the next section. 



It is worth pointing out that we did not assume any particular universality class for 
the transitions analyzed here in order to calculate their critical temperatures, since the 
analysis relied on the order parameter cumulant's properties. However, a systematic 



analysis of critical exponents was carried out in Ref. 18(, where the universality class was 
shown to be that of the 2D Ising model for square lattices and the two-dimensional Potts 
model with q = 3 for triangular lattices. 



III. ANALYTICAL APPROXIMATIONS AND COMPARISON BETWEEN SIMULATED 

AND THEORETICAL RESULTS 



Next, the transition is analyzed from the main theoretical models developed to treat 
the polymers adsorption problem. The study allows us to calculate the dependence of 
9 c {k) on k from three well-known multisite-adsorption theories: the Flory-Huggins (FH) 
approximation^, the Guggenheim-DiMarzio (GD) approach 4,5 and the recently devel- 
oped Semiempirical Model for Adsorption of Polyatomics (SE)~>2£. The corresponding 
expressions for the configurational entropy per site are: 
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From a given value of k (which depends on the approximation considered), the 2D 
[eqs. (EHUD] an d ID [eq. ([5])] curves cross at intermediate densities and two well differ- 
entiated regimes can be observed. In the first regime, which occurs at low densities, the 
2D approaches predict a larger entropy than the ID data. In the second regime (at high 
densities) the behavior is inverted and the 2D data present a smaller sjk^ than the ID 
results. Given that the theoretical results in 2D assume isotropy in the adlayer (interested 
readers are referred to Ref. 5], where this point is explicitly considered), the crossing of 
the curves shows that, in the second regime, the contribution to the 2D entropy from 
the isotropic configurations is lower than the contribution from the aligned states. Then, 
the existence of an intersection point is indicative of a I-N transition and allows us to 
estimate k m i n and 9 C from the different approximations studied. This intersection point 
can be easily calculated through a standard computing procedure; in our case, we used 
Mathematica software. The results are shown in Fig. 5 for the different studied cases. 
The symbology is indicated in the figure. 

In the case of square lattices, FH and GD predict values of k m i n = 3 and k min = 4, 
respectively. On the other hand, SE performs significantly better than the other ap- 
proaches, predicting the "exact" value of k min = 7™. With respect to triangular lattices, 
the values obtained for the critical k are k m i n = 4 (FH), k min = 4 (GD) and k m i n = 10 
(SE). Finally, the behavior of 9 C as a function of k for the analytical approximations fol- 
lows qualitatively the Monte Carlo simulation results, which reinforces the robustness of 
the analysis introduced here. 



IV. CONCLUSIONS 

In the present work, we have addressed the critical properties of long straight rigid 
rods on square and triangular lattices at intermediate density. The results were ob- 
tained by combining Monte Carlo simulations, finite-size scaling techniques and three well- 
known multisite-adsorption theories: the Flory-Huggins approximation, the Guggenheim- 
DiMarzio approach and the Semiempirical Model for Adsorption of Polyatomics. 

Two main conclusions can be drawn from the present work. On one hand, the critical 
density dependence on the particle size k has been reported. We found that 9 c (k) follows 
a power law as 9 c (k) oc k~ x . On the other hand, our analysis allowed us also to determine 
the minimum value of k {k m in = 7), which allows the formation of a nematic phase on a 
triangular lattice. 
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Figure Captions 

Fig. 1: Configurational entropy per site (in units of ks) versus surface coverage for 
different k as indicated and L/k = 20. Dotted line and symbols represent MC results 
for triangular lattices, and solid lines correspond to exact results for one-dimensional 
systems. Inset: Surface coverage dependence of the nematic order parameter for different 
k as indicated and L/k = 20. 

Fig. 2: (a) Curves of Ul vs 9 for square lattices, k = 10 and different lattice sizes as 
indicated. From their intersections one obtains 9 C . Inset: Size dependence of the order 
parameter, 5, as a function of coverage, (b) Same as part (a) for triangular lattices. 
Fig. 3: Simulated results for critical density 9 C as a function of k. The symbology is 
indicated in the figure. 

Fig. 4: (a) An example of a "small window" of side I = 8, which is embedded in an infinite 
lattice. Open squares correspond to empty sites and black bars represent adsorbed fc-mers 
with k = 4. (b) Same as part (a) for rescaled values of I and k. In this case, / = 24 and 
k = 12. 

Fig. 5: Theoretical results for critical density 9 C as a function of k. The symbology is 
indicated in the figure. 
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